## function

rm(list=ls())
library(ctgt)
library(stringr)
library(xlsx)

weed <- function(lams, accuracy = 50) {
  lams = lams[lams != 0]
  if (missing(accuracy)) {
    thresh <- 1/50
  } else {
    thresh <- 1/accuracy
  }
  lams <- -sort(-lams)
  qq <- length(lams)
  while ((qq>2 && lams[1] > 0) && (lams[qq] / lams[1] < thresh)) {
    q <- qq-1
    r <- qq-2
    lams[q] <- lams[q] + lams[qq]
    while ((r > 0) && (lams[r] < lams[q])) {
      lams[r:q] <- mean(lams[r:q]) 
      r <- r - 1
    }
    qq <- q
  }
  lams[1:qq]
}

cp = function(lam1, lam2){
  lam1= weed(lam1)
  lam2= weed(lam2)
  cdf1 = 0
  n1 = length(lam1)
  for(i in 1:n1){
    set.seed(123+i)
    cdf1 = cdf1+lam1[i]*rchisq(1e+6,1)
  }
  cdf2 = 0
  n2 = length(lam2)
  for(i in 1:n2){
    set.seed(1234+i)
    cdf2 = cdf2+lam2[i]*rchisq(1e+6,1)
  }
  cp1 = ecdf(cdf1)
  cp2 = ecdf(cdf2)
  cx = uniroot(function(x) cp1(x)-cp2(x), c((n1*0.05+1) * mean(lam2), (max(n1,n2)*0.05+1) * max(lam2)+max(lam2) ),extendInt = "yes" )$root
  alpha0 = 1-cp1(cx)
  return(alpha0)
  
}
##
##
##

load("d_uri.RData")
y = d_uri[,64]
X = d_uri[,1:63]

n = 77
m = 63


xs = str_replace_all(paste(rep("x",m),1:(m) ),fixed(" "), "") 
colnames(X) <- xs

sqrW = sqrt(mean(y)*(1-mean(y)) )# covariance of y
Z = X # full design matrix without intercept

WIHZ = sqrW *(sweep(Z,2,colMeans(Z))) ## w^{1/2}*(I-H)Z
IHZ = WIHZ/sqrW ##(I-H)Z



s1 = 1
set.seed(1234)
Lamf = round(eigen(WIHZ[,xs,drop=F]%*%t(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) #
hyps = str_replace_all(paste(rep("x",s1), sample.int(m, s1)),fixed(" "), "") 
Lams = round(eigen(WIHZ[,hyps,drop=F]%*%t(WIHZ[,hyps,drop=F]),symmetric = T,only.values = T)$values,8)  # ful


ca = c()
lv = c()
rxs = setdiff(xs,hyps)
for(rr in 1:(length(rxs)-1)){
  for(r in 1:8){
    set.seed((1+rr+r))
    HP = c(hyps,sample(rxs, rr) )
    LamHP = round(eigen(WIHZ[,HP,drop=F]%*%t(WIHZ[,HP,drop=F]),symmetric = T,only.values = T)$values,8)  
    hatlam = getL(Lamf,Lams,sum(LamHP))
    ca[(rr-1)*8+r] = cp(LamHP, hatlam)
    lv[(rr-1)*8+r] = sum(LamHP)
  }
}


write.xlsx(cbind(ca,lv),"ca1.xlsx")

####################################################################################
s1 = 10
set.seed(1234)
Lamf = round(eigen(WIHZ[,xs,drop=F]%*%t(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) #
hyps = str_replace_all(paste(rep("x",s1), sample.int(m, s1)),fixed(" "), "") 
Lams = round(eigen(WIHZ[,hyps,drop=F]%*%t(WIHZ[,hyps,drop=F]),symmetric = T,only.values = T)$values,8)  # ful


ca = c()
lv = c()
rxs = setdiff(xs,hyps)
for(rr in 1:(length(rxs)-1)){
  for(r in 1:10){
    set.seed((1+rr+r))
    HP = c(hyps,sample(rxs, rr) )
    LamHP = round(eigen(WIHZ[,HP,drop=F]%*%t(WIHZ[,HP,drop=F]),symmetric = T,only.values = T)$values,8)  
    hatlam = getL(Lamf,Lams,sum(LamHP))
    ca[(rr-1)*10+r] = cp(LamHP, hatlam)
    lv[(rr-1)*10+r] = sum(LamHP)
  }
}


write.xlsx(cbind(ca,lv),"ca10.xlsx")

####################################################################################
s1 = 30
set.seed(1234)
Lamf = round(eigen(WIHZ[,xs,drop=F]%*%t(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) #
hyps = str_replace_all(paste(rep("x",s1), sample.int(m, s1)),fixed(" "), "") 
Lams = round(eigen(WIHZ[,hyps,drop=F]%*%t(WIHZ[,hyps,drop=F]),symmetric = T,only.values = T)$values,8)  # ful


ca = c()
lv = c()
rxs = setdiff(xs,hyps)
for(rr in 1:(length(rxs)-1)){
  for(r in 1:16){
    set.seed((1+rr+r))
    HP = c(hyps,sample(rxs, rr) )
    LamHP = round(eigen(WIHZ[,HP,drop=F]%*%t(WIHZ[,HP,drop=F]),symmetric = T,only.values = T)$values,8)  
    hatlam = getL(Lamf,Lams,sum(LamHP))
    ca[(rr-1)*16+r] = cp(LamHP, hatlam)
    lv[(rr-1)*16+r] = sum(LamHP)
  }
}


write.xlsx(cbind(ca,lv),"ca30.xlsx")

####################################################################################
s1 = 50
set.seed(1234)
Lamf = round(eigen(WIHZ[,xs,drop=F]%*%t(WIHZ[,xs,drop=F]),symmetric = T,only.values = T)$values,8) #
hyps = str_replace_all(paste(rep("x",s1), sample.int(m, s1)),fixed(" "), "") 
Lams = round(eigen(WIHZ[,hyps,drop=F]%*%t(WIHZ[,hyps,drop=F]),symmetric = T,only.values = T)$values,8)  # ful


ca = c()
lv = c()
rxs = setdiff(xs,hyps)
for(rr in 1:(length(rxs)-1)){
  for(r in 1:30){
    set.seed((1+rr+r))
    HP = c(hyps,sample(rxs, rr) )
    LamHP = round(eigen(WIHZ[,HP,drop=F]%*%t(WIHZ[,HP,drop=F]),symmetric = T,only.values = T)$values,8)  
    hatlam = getL(Lamf,Lams,sum(LamHP))
    ca[(rr-1)*30+r] = cp(LamHP, hatlam)
    lv[(rr-1)*30+r] = sum(LamHP)
  }
}


write.xlsx(cbind(ca,lv),"ca50.xlsx")
